#实践2####

#建立中文对象，英文对象，隐藏对象
English="CSU"
中文="中南大学"
.hidden="生命科学学院"
#建立数值向量和因子，利用tapply函数计算均值
vector<-c(1:10)
factor <- factor(rep(c("Group1", "Group2", "Group3"), times = c(3, 3, 4)))
means <- tapply(vector,factor,mean)
#建立数值向量，与上述向量和因子合并成数据框
vector1<-(11:20)
vector2<-(21:30)
data<-data.frame(Num = vector, Num1 = vector1, Num2 = vector2,Group = group_factor)

#实践3####
#实践3第一题

#ADdata数据读取
data1<-read.table("D:/c dongzhiying/MEDownloads/ADdata/ADdata1.txt")
data2<-read.csv("D:/c dongzhiying/MEDownloads/ADdata/ADdata2.csv")
data2<-read.csv("D:/c dongzhiying/MEDownloads/ADdata/ADdata2.csv",row.names = 1)
library("openxlsx")
data3<-read.xlsx("D:/c dongzhiying/MEDownloads/ADdata/ADdata3.xlsx",rowNames = TRUE)
data4<-read.table("D:/c dongzhiying/MEDownloads/ADdata/ADdata4.txt")

#对列名蛋白质ID排序
common_ids <- Reduce(intersect, list(rownames(data1), rownames(data2),rownames(data3),rownames(data4)))
data1 <- data1[common_ids, ]
data2 <- data2[common_ids, ]
data3 <- data3[common_ids, ]
data4 <- data4[common_ids, ]

#合并数据
merge_data1 <- cbind(data1,data2)
merge_data2 <- cbind(data3,data4)
merge_data3 <- cbind(merge_data1,merge_data2)

#保存为csv,xlsx,txt格式文件
write.csv(merge_data3, "merged_data.csv" )
write.xlsx(merge_data3, "merged_data.xlsx" )
write.table(merge_data3, "merged_data.txt")

#实践3第二题

#读取文件
clin_inf<-read.csv("D:/c dongzhiying/MEDownloads/class2/clin_inf.csv",sep =" ",)
count<-read.csv("D:/c dongzhiying/MEDownloads/class2/count.csv",sep=" ")
exp_inf<-read.csv("D:/c dongzhiying/MEDownloads/class2/exp_inf.csv",sep=" ")
ID_annoation<ad.csv("D:/c dongzhiying/MEDownloads/class2/ID_annoation.csv",sep=" ")

#数据整理

#合并临床数据和实验数据
meta <- merge(clin_inf,exp_inf)
#修改行名
rownames(count) <- count$gene_id 
count <- count[,-1] 
#加载所需包
install.packages("BiocManager")
BiocManager::install("Biobase")
library(BiocGenerics)
library(Biobase)
#规范化整理样本，修改列名
x <- colnames(count)
x <- gsub("(.*).$","\\1",x)
x <- gsub("X(...$)","PTB\\1",x)
x <- gsub("X","XYA",x)
colnames(count) <- x
#数据对齐
rownames(meta) <- meta$样本名称
meta <- meta[match(colnames(count),rownames(meta)),]
#简化基因ID，去除重复基因并对齐表达矩阵
rownames(ID_annoation) <- ID_annoation$gene_id
ID_annoation$gene_id <- gsub("\\..*","",ID_annoation$gene_id)
pos2 <- duplicated(ID_annoation$gene_id)
ID_annoation2 <- ID_annoation[ID_annoation$gene_id %in% ID_annoation$gene_id[pos2], ] 
ID_annoation <- ID_annoation[!duplicated(ID_annoation$gene_id), ]
rownames(ID_annoation) <- ID_annoation$gene_id
ID_annoation <- ID_annoation[match(rownames(count), rownames(ID_annoation)), ]
#构建ExpressionSet对象
eset <- ExpressionSet(
  assayData = as.matrix(count),
  phenoData = new("AnnotatedDataFrame", data = meta),
  featureData = new("AnnotatedDataFrame", data = ID_annoation)
)
#提取/删选部分基因
eset2 <- eset[1:10, 1:10]
treated_samples <- eset[, eset$批次 == "batch1"]
chrX_genes <- eset[fData(eset)$seqnames == "chr1", ]

#实践3第三题

num=sample(1:50,10,replace = TRUE)

#for循环
count <- 0
for (i in num) {
  if (i > 10) {
    count <- count + 1
  }
}
print(paste("大于10的元素个数:", count))

#while循环
count <- 0
index <- 1
while (index <= length(num)) {
  if (num[index] > 10) {
    count <- count + 1
  }
  index <- index + 1
}
print(paste("大于10的元素个数:", count))



#实践4_画图####

plot(c(5,4,3,20,1))
plot(2,10)
plot(-2,10)

pdf("1.pdf")
#本处写绘图代码
dev.off()
#箱图
boxplot(merge_data3[1:10,1:10])
#条形图
barplot(as.matrix(merge_data3[1:5,1:5]))
x<-merge_data3[1:5,1:5]
#基于基因绘制直方图

rownames(count) <- count$gene_id 
count <- count[,-1] 
hist(count[,1])
hist(log10(count[,1]))
hist(log10(count[,1]+1))

hist(count[1,])
a<-as.matrix(count[1,])
a<-unlist(count[1,])
hist(a,breaks = 20)
pos=which(genemed>100&genemed)
hist(unlist(count[9,]),breaks = 20)
hist(unlist(count[pos[1],]),breaks = 20)


hist(count$PTB008D,breaks = 100)
hist(log2(count$PTB008D+1))
hist(log2(count$PTB008D+1),breaks = 100)
hist(log2(count$PTB080D+1),breaks = 50)
#基于样本绘制直方图
hist(log(merge_data3$asym54+1),breaks=50)
hist(log(merge_data3$asym26+1),breaks=50)
#图形美化
hist(log2(count$PTB008D + 1),
     breaks = 114,
     main = "PTB008D样本表达量分布",
     xlab = "Log2(表达量 + 1)", 
     ylab = "基因频率",
     col = "#69b3a2",
     border = "white",
     cex.lab = 1,
     cex.main = 1)

#pca分析
pca <- prcomp(t(merge_data3))
x <- data.frame(t(merge_data3))
plot(x$O00451,x$O00161)

install.packages("ggrepel")
library(ggplot2)
library(ggrepel)
x1 <- data.frame(pca$x)
x1$label <- rownames(x1)
ggplot(x1,aes(PC1,PC2)) +
  geom_point() +
  geom_label_repel(aes(label = label))
plot(pca$x[,1],pca$x[,2])
plot(pca)
biplot(pca,var.axes=FALSE)

#聚类分析
d = dist(t(merge_data3))
clu <- hclust(d)
plot(clu)

#火山图
load("D:/c dongzhiying/MEDownloads/volcano.RData")
#数据处理
plotdata <- prostat
#列名
colnames(plotdata)[2] <- c("log2fc")
plotdata$group <- "nosig"
#义差异表达分组
pos = plotdata$log2fc>0.58 & plotdata$P<0.05
plotdata$group[pos] <- "up"
pos = plotdata$log2fc< -0.58 & plotdata$P<0.05
plotdata$group[pos] <- "down"
plotdata$label <- plotdata$ID
plotdata$label[plotdata$group == "nosig"] <- ""
library(ggplot2)
#保存火山图.jpg
jpeg("volcano.jpg")
#绘图
ggplot(plotdata,aes(x=log2fc,y=-log10(P))) +
  geom_point(aes(color = group)) +
  scale_color_manual(values = c("green","grey","red"),
                     limits = c("down","nosig","up")) +
  geom_label_repel(aes(label = label))
dev.off()

jpeg("volcano1.jpg")
x <- plotdata;
x$log10p <- -log10(x$P)
x$color <- x$group
x$color[x$color == "nosig"] <- "grey"
x$color[x$color == "up"] <- "red"
x$color[x$color == "down"] <- "green"
#图雏样
plot(x$log2fc, x$log10p, type = "n",
     main="volcano plot",xlab="FC",ylab="-log10 p-value")
#加点
points(x$log2fc,x$log10p,col=x$color,pch=16,cex=1)
points(1.5,6,col="black")
#加线
abline(v=log2(1.5),lty=3)
abline(v=log2(1/1.5),lty=3)
abline(h=-log10(0.05),lty=3)
#加图例
legend(-1.5,1,c("up","nosig","down"),col = c("red","grey","green"),text.col ="black",
       pch = c(16,16,16),cex = c(0.4,0.4,0.4))

#在图中标出up基因
up <- subset(x,color == "red")
text(up$log2fc,up$log10p+0.2,labels = up$ID)
dev.off()



#实践5 R统计分析####
#方差分析#
ADdata<-merge_data3
# 数据对数化处理
pos = ADdata == 0;
ADdata[pos] <- NA;
ADdata[1,1]
ADdata2 <- log2(ADdata)
colnames(ADdata)
# 提取分组信息
group <- gsub("\\d","",colnames(ADdata))
table(group)

# 提取单个蛋白数据测试
pro1 <- ADdata2[i,];
pro1.1 <- unlist(pro1);
pro1.2 <-as.vector(t(pro1))
pro1.1
# 判断各组有效数据量
# 提取各组数据
posctl <- group == "ctl"
ctl <- pro1.1[posctl]
posad <- group == "ad"
ad <- pro1.1[posad]
posasym <- group == "asym"
asym <- pro1.1[posasym]
# 计算各组非NA的个数
posx <- !is.na(ctl)
ctlnum = sum(posx)
posx2 <- !is.na(ad)
adnum = sum(posx2)
posx3 <- !is.na(asym)
asymnum = sum(posx3)
# 方差分析或输出NA
if(asymnum >3 & adnum>3 & ctlnum >3){
  anovaresult <- oneway.test(pro1.1~group)
  p <- anovaresult$p.value
} else p = NA
# 初始化结果存储向量
protein_names <- rownames(ADdata2)
p_values <- vector("numeric", length = length(protein_names))
names(p_values) <- protein_names
# 代入循环
# 循环处理每个蛋白
for (i in 1:nrow(ADdata2)) {
  # 提取当前蛋白数据
  pro1 <- ADdata2[i, ]
  pro1.1 <- unlist(pro1)
  
  # 提取各组数据
  posctl <- group == "ctl"
  ctl <- pro1.1[posctl]
  
  posad <- group == "ad"
  ad <- pro1.1[posad]
  
  posasym <- group == "asym"
  asym <- pro1.1[posasym]
  
  # 计算各组非NA的个数
  ctlnum <- sum(!is.na(ctl))
  adnum <- sum(!is.na(ad))
  asymnum <- sum(!is.na(asym))
  
  # 进行方差分析或赋值为NA
  if (asymnum >= 3 & adnum >= 3 & ctlnum >= 3) {
    anovaresult <- oneway.test(pro1.1 ~ group)
    p_values[i] <- anovaresult$p.value
  } else {
    p_values[i] <- NA
  }
}
# 整理结果
result_df <- data.frame(
  Protein = protein_names,
  P_value = p_values,
  stringsAsFactors = FALSE
)
# 添加多重检验校正
result_df$Adj_P_value <- p.adjust(result_df$P_value, method = "BH")
# 查看结果
head(result_df)
summary(result_df$P_value)
#P值分布直方图
# 过滤掉NA值
valid_p <- na.omit(result_df$P_value)
# 绘制P值分布
hist(valid_p, breaks = 30, col = "skyblue", 
     main = "Distribution of ANOVA P-values",
     xlab = "P-value", ylab = "Frequency")
abline(v = 0.05, col = "red", lty = 2)
legend("topright", legend = "P = 0.05", col = "red", lty = 2)
#富集分析###
# 加载包
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("enrichplot")
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)

install.packages("gprofiler2")
library(gprofiler2)
library(ggplot2)

install.packages("enrichR")
library(enrichR)
#load("D:/MyDownload/Rstudio/dzy_rwork/vovolcano.RData")
load("D:/c dongzhiying/MEDownloads/volcano.RData")
# 提取显著差异表达蛋白P < 0.05
sig_genes <- prostat[prostat$P < 0.05 & !is.na(prostat$P), "ID"]
print(paste("显著基因数量:", length(sig_genes)))
sig_proteins <- subset(volcano, P < 0.05)$ID
# 选择数据库（GO Biological Process 2023）
dbs <- c("GO_Biological_Process_2023")
# 执行富集分析
enrich_results <- enrichr(sig_genes, dbs)
# 提取结果
go_results <- enrich_results[[1]]  # 第一个数据库的结果
# 筛选显著富集项（调整后P < 0.05）
sig_go <- subset(go_results, P.value < 0.05)
# 按P值排序并保留前20项
top_go <- head(sig_go[order(sig_go$P.value), ], 20)
# 简化Term名称（可选）
top_go$Term <- gsub(" \\(GO:\\d+\\)", "", top_go$Term)
library(ggplot2)
ggplot(top_go, aes(x = -log10(Adjusted.P.value), 
                   y = reorder(Term, -log10(P.value)),
                   size = -log10(P.value),
                   color = Odds.Ratio)) +
  geom_point(alpha = 0.7) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "GO Biological Process 富集分析",
       x = "-Log10(Adj.P)", 
       y = "",
       color = "Odds Ratio",
       size = "-Log10(Adj.P)") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 10))
# 加载包
BiocManager::install("Biostrings")
BiocManager::install("DOSE")
BiocManager::install("GO.db")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("enrichplot")
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
# 提取显著差异表达蛋白（假设数据框包含Gene、P.Value和logFC列）
sig_genes <- prostat[prostat$P.Value < 0.05 & !is.na(prostat$P.Value), "ID"]
# 提取显著差异表达蛋白P < 0.05
sig_genes <- prostat[prostat$P < 0.05 & !is.na(prostat$P), "ID"]
print(paste("显著基因数量:", length(sig_genes)))
#sig_proteins <- subset(volcano, P < 0.05)$ID
# 将基因符号转换为Entrez ID（假设基因名为symbol）
gene.df <- bitr(sig_genes, fromType = "SYMBOL",
                toType = c("ENTREZID"),
                OrgDb = org.Hs.eg.db)

entrez_ids <- gene.df$ENTREZID
# GO富集分析（生物学过程）
ego <- enrichGO(gene          = entrez_ids,
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",  # 生物学过程
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.2,
                readable      = TRUE)

# 查看富集结果
head(ego)
# 简单气泡图
dotplot(ego, showCategory=15) + 
  ggtitle("GO Biological Process Enrichment") +
  theme(plot.title = element_text(hjust = 0.5))

# 升级气泡图
jepg("go_enrichment.jpg")
png("GO_Enrichment.png")
ggplot(ego@result[1:15, ], aes(x=GeneRatio, y=Description, 
                               size=Count, color=-log10(p.adjust))) +
  geom_point(alpha=0.8) +
  scale_color_gradient(low="blue", high="red") +
  labs(title="GO Enrichment Analysis (Top 15 Terms)",
       x="Gene Ratio", 
       y="GO Term",
       color="-log10(Adj.P)",
       size="Gene Count") +
  theme_minimal() +
  theme(axis.text.y = element_text(size=10),
        plot.title = element_text(hjust=0.5, face="bold"))
dev.off()


#实践6 转录组分析####
#加载包
library(BiocGenerics)
library(Biobase)
library(matrixStats)
library(MatrixGenerics)
library(GenomicRanges)
library(stats4)
library(S4Vectors)
library(IRanges)
library(SummarizedExperiment)
library(GenomeInfoDb)
BiocManager::install("DESeq2")
library(DESeq2)
install.packages("tidyverse")
library(tidyverse)

library("DESeq2")
BiocManager::install(c("limma", "edgeR"))
library("limma")
library("edgeR")
BiocManager::install(c( "pheatmap","RColorBrewer"))
install.packages("biomaRt")
BiocManager::install("biomaRt")
library(biomaRt)
library(pheatmap)
library(RColorBrewer)
#四舍五入整数化
count1 <- round(count)
#匹配样本
shared_samples <- intersect(colnames(count1), clin_inf$样本名称)
count1 <- count1[, shared_samples]
clin_inf1 <- clin_inf[clin_inf$样本名称 %in% shared_samples, ]
# 确保样本顺序一致
clin_inf1 <- clin_inf1[match(shared_samples, clin_inf1$样本名称), ]
rownames(clin_inf1) <- clin_inf1$样本名称
#分组
clin_inf1$age_group <- ifelse(clin_inf1$年龄 < 65, "middle", "old")
clin_inf1$age_group <- factor(clin_inf1$age_group, levels = c("middle", "old"))
# 统计各组样本数
table(clin_inf1$age_group)
# 创建DESeqDataSet对象
dds <- DESeqDataSetFromMatrix(
  countData = count1,          # 原始计数矩阵
  colData = clin_inf1,             # 临床数据（含age_group）
  design = ~ age_group             # 设计公式
)

# 过滤低表达基因（至少在10%样本中count≥10）
keep <- rowSums(counts(dds) >= 10) >= ceiling(0.1 * ncol(dds))
dds <- dds[keep, ]

# 检查样本分组分布
table(colData(dds)$age_group)

# 可视化分组情况
library(ggplot2)
ggplot(as.data.frame(colData(dds)), aes(x=age_group, fill=age_group)) +
  geom_bar() +
  labs(title="Sample Distribution by Age Group")
# 标准分析流程
dds <- DESeq(dds)
# 提取结果（老年组 vs 中年组）
res <- results(dds, 
               contrast = c("age_group", "old", "middle"),
               alpha = 0.05,        # 显著性阈值
               lfcThreshold = 0)    # 不设log2FC阈值
# 添加基因符号
library(AnnotationDbi)
library(org.Hs.eg.db)
res$symbol <- mapIds(org.Hs.eg.db,
                     keys = rownames(res),
                     column = "SYMBOL",
                     keytype = "ENSEMBL")
# 按调整p值排序
res_ordered <- res[order(res$padj), ]
# 筛选显著差异基因（padj < 0.05 & |log2FC| > 1）
sig_genes <- subset(res_ordered, padj < 0.05 & abs(log2FoldChange) > 1)
# 保存结果
write.csv(as.data.frame(res_ordered), "all_genes_age_diff.csv")
write.csv(as.data.frame(sig_genes), "significant_genes_age_diff.csv")
# 查看Top差异基因
head(sig_genes, 10)

# 火山图
install.packages("EnhancedVolcano")
BiocManager::install("EnhancedVolcano")
library(ggrepel)
library("EnhancedVolcano")
EnhancedVolcano(res,
                lab = res$symbol,
                x = 'log2FoldChange',
                y = 'padj',
                pCutoff = 0.05,
                FCcutoff = 1,
                title = "Old vs Middle Age (≥65 vs <65)")

png("volcano_plot.png", width = 800, height = 600, res = 100)
ggplot(as.data.frame(res_ordered), aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(alpha = 0.6, size = 2, 
             aes(color = ifelse(padj < 0.05 & abs(log2FoldChange) > 1, 
                                ifelse(log2FoldChange > 1, "Up-regulated", "Down-regulated"), 
                                "Not significant"))) +
  scale_color_manual(values = c("Not significant" = "gray", 
                                "Up-regulated" = "red", 
                                "Down-regulated" = "blue")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray50") +
  labs(title = "火山图：老年组 vs 中年组差异表达基因",
       x = "Log2 Fold Change",
       y = "-Log10(Adjusted p-value)",
       color = "基因类别") +
  theme_minimal() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 10))
dev.off()
#cat("差异表达基因统计：\n")
#cat("上调基因 (padj < 0.05 & log2FC > 1):", nrow(subset(sig_genes, log2FoldChange > 1)), "\n")
#cat("下调基因 (padj < 0.05 & log2FC < -1):", nrow(subset(sig_genes, log2FoldChange < -1)), "\n")
#cat("总计显著差异基因:", nrow(sig_genes), "\n") 
